Density-operator theory of orbital magnetic susceptibility in periodic insulators 
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The theoretical treatment of homogeneous static magnetic fields in periodic systems is chaUenging, 
as the corresponding vector potential breaks the translational invariance of the Hamiltonian. Based 
on density operators and perturbation theory, we propose, for insulators, a periodic framework 
for the treatment of magnetic fields up to arbitrary order of perturbation, similar to widely used 
schemes for electric fields. The second-order term delivers a new, remarkably simple, formulation 
of the macroscopic orbital magnetic susceptibility for periodic insulators. We validate the latter 
expression using a tight-binding model, analytically from the present theory and numerically from 
the large-size limit of a finite cluster, with excellent numerical agreement. 

PACS numbers; 71.15.-m,75.20.-g 
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I. INTRODUCTION 

The ability to compute the response of periodic sys- 
tems to homogeneous electric fields, strain and atomic 
displacements is a key ingredient in our current under- 
standing of dielectric materials (also ferroelectrics, piezo- 
electrics): from first principles, '^"'^ one obtains easily 
the polarization, dielectric constants, piezoelectric coef- 
ficients, phonon band structure, etc. Compared to the 
treatment of such responses for molecules, the handling of 
periodic boundary conditions has raised numerous chal- 
lenges: for instance, the linear potential associated with 
a static homogeneous electric field breaks the transla- 
tional symmetry. While such problems have been suc- 
cessfully addressed for the above-mentioned perturba- 
tions, the treatment of homogeneous magnetic fields is 
not as mature. With the current interest in multiferroi'c 
materials)^ and the long-term interest in magnetic field- 
based spectroscopies, a unified framework for all these 
responses is highly desirable. 

The first strategy followed to treat a periodicity break- 
ing was to consider perturbations with specific com- 
mensurate wavevectors and corresponding supercells. 
For atomic displacements, this approach is known as 
the "frozen-phonon" method.^ Homogeneous magnetic 
fields can be treated in this spirit, with an artificial 
modulation," although (1) working with supercells is 
CPU time-consuming, (2) the study of couplings is te- 
dious in this approach. More powerful formalisms, based 
on perturbation theory, that do not rely on supercells or 
long- wavelength limits, have been developed for atomic 
displacements and electric fields j^i^ii^""— For phonons, the 
long- wavelength phases can be factorized, such that a 
purely periodic treatment is recovered. For the electric 
field, the position operator can be replaced by the dif- 
ferentiation with respect to the wavevector. The Berry 
phase approach to the electrical polarization is probably 
the most striking consequence of the latter link.— Owing 
to these advances, linear and nonlinear responses can be 



addressed, as well as couplings between different pertur- 
bations, in a purely periodic primitive cell framework. 

For homogeneous magnetic fields in insulators the diffi- 
culties outlined above are more severe. We focus only on 
the orbital coupling; coupling to spin is not affected by 
periodicity issues. The presence of the magnetic field not 
only breaks the periodicity of the Hamiltonian but also 
induces a vector coupling to the electron dynamics. In a 
pioneering work, Mauri and Louie proposed a method to 
compute the orbital magnetic susceptibility (OMS) from 
the long wavelength limit of an oscillating perturbationj^ 
The theory has been adapted to the computation of the 
chemical shielding tensor and related quantities.—"— Al- 
though relying on perturbation theory concepts to avoid 
the use of supercells, this formalism introduces auxiliary 
oscillating quantities that are not consistent with the pe- 
riodicity of the lattice and that break the rotational in- 
variance of the global system. Thus for example the ten- 
sorial structure of the expression of the OMS could not be 
recovered.^ In a related alternative approach, by Sebas- 
tiani and Parrinello,— the position operator is replaced 
by localized sawtooth potentials, one for each orbital. In 
practice, the use of supercells, needed to deal with the 
spatially localized Wannier functions, is still necessary in 
this approach. 

Recently, a theory of orbital magnetization has been 
proposed, in which only periodic Bloch wavefunctions 
and Hamiltonian are used,^** bringing the understanding 
of this (bulk) quantity to the same formal level as the 
electric polarization. Based on this result, the orbital 
magnetoelectric coupling has been derived from density- 
matrix perturbation theoryJ^ No use of a supercell or 
long-wavelength limit is needed in this approach. 

In the present contribution, we show how to apply the 
density operator approach of Ref. [isl . to arbitrary or- 
ders of perturbation in the magnetic field. Focusing on 
the second order in this expansion, we obtain a new for- 
mula for the OMS, based on the first-order response of 
the density operator to magnetic field and wavevector. 
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in a purely periodic framework. The present approach is 
compatible with the similar treatment of electric fields, of 
atomic displacements and of their couplings, to arbitrary 
orders. We then check the theory by considering a two- 
dimensional (2D) periodic tight-binding (TB) model, for 
which, thanks to the new approach, the OMS can be ob- 
tained as an integral over the Brillouin Zone (BZ) of an 
analytical expression. Alternatively, we numerically solve 
this model for clusters of increasing size, considering ex- 
plicitly the magnetic field. Essentially exact agreement 
is obtained between the two approaches. 

Alternative OMS formulae were proposed already fifty 
years ago,^^ in the context of effective Hamiltonians or 
TB models. However, none of them seem compatible 
with the currently used formalisms for electric or atomic 
displacement responses. The present approach is actually 
considerably simpler. Furthermore, to our knowledge, no 
OMS formula for a periodic insulator has ever been vali- 
dated by comparison with numerical results on a solvable 
cluster model in the large size limit. 



II. THEORY 



Let us recall the first steps in the approach of Ref. |15|. 
In atomic units, the Hamiltonian H for an electron in 
a vector potential A(r) is H — ^ zV - ^A) -I- V{r), 
where c is the speed of light, and the potential V is pe- 
riodic for lattice vectors R, that is, V{r + R) = V{r). 
Choosing the gauge A = x r, the Hamiltonian does 
not possess the translation symmetry of the lattice, but 
does have magnetic translation symmetryJ^ The kernel 
of operators O possessing this symmetry can be related 
to a periodic kernel O: 



where 



0(ri,r2) = (5(ri,r2)e 



-iB-ri Xr2/2c 



(1) 



with 0(ri -I- R,r2 -I- R) = 0(ri,r2). Crucially, the pe- 
riodic counterpart of the Hamiltonian in this approach 
has no vector potential dependence (nor magnetic field 
dependence): H = -|V^ + V{r). 

The density operator corresponding to an Hamilto- 
nian, for an insulator at Kelvin, can be obtained by 
the minimization of the expectation value of the energy 
of the system on the ensemble of idempotent density ma- 
trices {i.e. p = pp) with fixed electron number: 



min {Tr[pi7(B)]}. 

p idempotent 



(2) 



At the minimum, [p, H] = 0. The density operator also 
possesses magnetic translation symmetry, so that Eq. ([2]) 
can be recast in terms of H and p provided that the 
product of the two operators is transformed according 
to Eq. dl]). In real space the product of two operators 
r =VW becomesi^iil 



f(ri,r3)= / dr2V(ri,r2)W(r2,r3)e-*"*i= 



(3) 



^123/ 



B.(ri X r2 + r2 X r3 +r3 X ri)/2 (4) 



is proportional to the magnetic fiux through triangle 123. 
Based on Eq. ([3]), p is no longer idempotent, due to the 
appearance of the e"*"^!^^/"*" phase factor. In Ref. [HI, p 
is expanded to first order in B, before considering the 
decomposition of the operators in the BZ, in order to 
obtain the orbital magnetoelectric coupling. 

We find that the idempotency problem can be avoided 
by transforming Eq. (|3]) to a combined BZ and primitive 
cell integral, before any expansion in B. To obtain this 
result, we decompose periodic operators 0(ri,r2) into 
operators that are separately periodic in each argument 
and characterized by a wavevector: 



0(ri,r2) 



BZ 



dk 

(2^ 



■e"^"-iOk(ri,r2)e-*''-^% 



(5) 



with Ok(ri +R,r2) = Ok (ri , r2 + R) = Ok(ri,r2), and 
Ck(ri,i"2) — C'k+G(ri,r2) for all reciprocal-lattice vec- 
tors G. Under this decomposition, Eq. ([3]) becomes 



r(ri,r3) = 



dk 



BZ (2^)3 jno 
xVk(ri,r2)yVk+Ak(B)(r2,r3) 



(6) 



where fig is the primitive cell volume and Ak — (r^ — 
ri) X B/2c is a position- and magnetic field- dependent 
increment of the wavevector k. Such an increment defines 
an ri-dependent kernel, and hence, an ri-dependent op- 
erator. Although unusual, these expressions are math- 
ematically well-defined. The transformation proceeds 
with the Taylor expansion of W with respect to B, fol- 
lowed by multiple integrations by parts with respect to 
k. This procedure yields the full expansion of Tk to all 
orders in B. For simplicity, we denote the derivatives 
with respect to k in direction a by omit real-space 
arguments, and use the summation convention with e the 
totally antisymmetric unit tensor. The expansion is then 



m / m 



m— 1 ^ ^ \n— 1 / 

xi^p^■■■^pMi^^,■■■^^^Wk),{7) 
and explicitly to second order, 
Tk = VkWk + ii/2c)ie^p^Bo,){dfjVk){d^Wk) 

(8) 

The density-matrix perturbation theory valid for arbi- 
trary orders of perturbation, as developed in k-space for 
homogeneous electric fields,— can now be generalized to 
magnetic fields. Defining perturbation orders by A such 
that B(A) = AB^^i, we expand a generic quantity X 
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as X{B) = + AX(i) + A2x(2) + . . .. The periodic 
counterpart of the Hamihonian, i/^ has no magnetic field 
dependence, as outhned before, hence, Hk = while 

for TO 7^ 0, =0. The density operator at any order 
can be decomposed into different blocks, acting inside 
the occupied subspace (denoted V for Valence), or inside 
the unoccupied subspace (C for Conduction), or coupling 
different subspaces CV/VC (we refer to the subspaces 
obtained at zero B). As Hk — H^'^ is block-diagonal in 
these subspaces, the expansion of the energy yields 



BZ 



dk 



Tr[(pL"^v+PSc)^k] 



(9) 



The first order term is related to the magnetization M, 
iJ*^^) = B • M, while the second order term is related to 
the QMS X, E(^'> =2^,^. B,x^JB,. 

The pI^"'' can be found recursively from the p'^^ with 
I < n, by an adaptation of the density-matrix pertur- 
bation theory of Ref. [H. The block-diagonal parts of 
density matrices are given by 



-in) 



~(0) -(n) ~(0) 



pL"^ = (i-pL"^)pL"z](i-Pr) 

where, for first and second orders, 

Pk£)= ^£^07^" (9/3/5^°^ )(a^p|^°^). 



(10) 

(11) 

(12) 



~(2) ~(1)~(1) 
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For the full p^^\ the off-diagonal CV blocks are needed 
in addition to the first-order CC and VV blocks, and are 
obtained by solving the equation 

x[(a^pf )(a^i?k) - idpH^)idJ^Wu\ (14) 

where the r.h.s expression is a CV projection of the 
angular-momentum operator times the magnetic field, in 
reciprocal spaceJ^ 

For n = 1, these equations reduce to Eqs. (21), (23), 
and (28) of Ref. [l^ where they have been further ex- 
pressed in terms of Bloch wavefunctions for occupied 
states and their first fc-derivatives. When substituted 
into the first-order term of Eq. ([9]), the multi-band for- 
mula derived in Ref. for the magnetization is obtained 
(as B • M). Likewise, the second-order equations can be 
expressed in terms of the Bloch wavefunctions for occu- 
pied states only, their first-order changes, and their first- 
and second- fc-derivatives, yielding an explicitly periodic 
formulation of the OMS. These rather lengthy expres- 
sions will be detailed elsewhere. The tensorial structure 



FIG. 1. Tight-binding model used to check the theory. A 
square lattice of A sites and B sites are coupled by energies t 
(A-B nearest-neighbors) and s (A- A nearest neighbors). 

A I 



of Eqs. (|9HT4|) is obvious, as the terms depend on n fac- 
tors of B in different directions. 

The OMS is seen from Eq. ([9]) to arise from both the 
valence and the conduction subspaces, like the magne- 
tization. For both subspaces, from Eq. (|13p . there are 
three contributions: a term quadratic in the density- 
operator response p'^\ a term linear in this response, and 
a term that is independent of the response of the elec- 
trons. The quadratic and linear CV contributions can be 
linked to each other on the basis of Eq. ([T4| . Their sum, 
being always negative and due to density matrix relax- 
ation, yields the Van Vleck paramagnetic contribution 
in the present formalism. Alternative decompositions 
of the OMS Eq. ([T^ exist, similarly to the different ex- 
pressions for the dielectric susceptibility in density func- 
tional perturbation theory, see, e.g., Eqs. (37) and (38) of 
Ref. i. The expression for the OMS presented here can 
be shown to be variational, and delivers Eq. (fT4|) from 
Euler-Lagrange conditions. 



III. VALIDATION OF THEORY 

The predictions of the theory were checked using a 
tight-binding model (Fig.[T]). This model consists of a 2D 
square lattice with A sites at the vertices and B sites at 
the center of each square, with on-site energies Ea < Eb. 
Nearest-neighbor A, B pairs are coupled by energy t, and 
nearest- neighbor A, A pairs by energy s (Fig. [J). The 
A sites are defined as initially occupied and the B sites 
empty. The two different types of couplings are necessary 
to exhibit non-trivial behavior from all three contribu- 
tions to the susceptibility; with only the more obvious t 
couplings, only the frozen-density-operator contribution 
is non-zero. Following Ref. [l^ the presence of a mag- 
netic field perpendicular to the plane of the model is in- 
cluded by multiplying the Hamiltonian matrix elements 
by e-^B(raXrb)/2 ^ where Tq and describe the locations 
of the coupled sites. 

This model is simple enough for Eqs. (|9HT4|) to lead 
to analytical expressions, apart for a global integral over 
the 2D BZ. The occupied (valence) band eigenenergies 
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-Bki) are given by 



with 



Ak 



— s(cos kx + cos ky), 
-At cos{kx/2) cos(fc,y/2). 



(15) 

(16) 

(17) 
(18) 



The periodic part of the Bloch eigenfunctions for the va- 
lence band are 



where 



\ukv) = cos6'k|A) + sin6'k|-B), 



tan0k = (l-^k-i?k)/Ak. 



(19) 



(20) 



With these definitions, the frozen, Unear, and quadratic 
density operator contributions to x are respectively 



Xfr, 



Xiii 



BZ 



4(27r)= 



dki dki 



-2Xquad = -2(27r) 



BZ 



(21) 



with 



i?k = .st(l + ^k)sin(^)sin(^) 
x[cos^(|)-cos^(|)]. 



(22) 



These integrals were carried out numerically to 10"^'' 
convergence. 

The susceptibility of the model was also computed di- 
rectly, by diagonalization of the Hamiltonian matrix for 
finite size lattices. For an iV x grid the energy con- 
verged roughly like 1/A^. For a given parameter triple 
(i, s, B) the energy was computed for a range of lat- 
tice sizes [N = 10- • -100), and the resulting values fit 
to a fourth order polynomial in 1/A^. This procedure 
yielded an estimate of the infinite-size limit accurate 
to about 1 part in 10^. To compute the second order 
energy change with magnetic field, energies as a func- 
tion of mesh size were computed for parameter triples 
(s,i, i?), where B was one of 0.00, 0.05, and 0.10, and 
s and t were fixed. Then was estimated from a 

finite-difference formula, valid when E{B) — E{—B): 
« [16£;(B/2)-£;(B)-15£;(0)]/(3B2), for B = 0.10. 

Fig. [5] shows the agreement between the theoretical 
prediction of E'^'^^ and that obtained by exact diagonal- 
ization. Note that at large s or i?*^^) has substantial 
contributions from the linear and quadratic components, 
in addition to the frozen term. In all cases agreement 
between theory and numerical diagonalization is on the 
order of a few parts in 10* or better. 

Beyond the OMS, the knowledge of the first-order den- 
sity matrix allows computation of every coupling between 



FIG. 2. Magnetic susceptibility for the tight-binding model 
with (a) t — 2.0 and a range of s values, and (b) s — 0.2 and 
a range of t values. In both cases, the full theory (solid line) 
is composed of three terms (dashed line: frozen wavefunction; 
dotted: linear in p'-^^; dash-dot: quadratic in p'^'). Squares: 
exact diagonalization values. 




a magnetic field and other (regular) perturbations of the 
system, as well as the orbital current as needed for the 
nuclear magnetic resonance (NMR) shielding. Indeed, 
the mixed derivative of the energy with respect to the 
magnetic field (indexed by A) and another perturbation 
(indexed by fi) is given by 



dk 



BZ 



(27r) 



(23) 



where, crucially, the derivative of the first-order density 
matrix with respect to /x is not needed. 



IV. CONCLUSIONS 

The present approach to the OMS through expansion 
of the total energy, Eq. ([9]), places it in a new frame- 
work, in which its link with the orbital magnetization 
(from the first-order term) is clear, and higher-order sus- 
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ceptibilities can be computed as well. Of course, the in- 
terest in higher-order derivatives of the total energy with 
respect to the magnetic field is purely academic. How- 
ever, electric fields and vibrational effects can be treated 
on the same footing, at linear and non-linear orders, so 
that the path is opened to a unified approach to coupled 
magnetic, electric and thermodynamic effects in insula- 
tors, expected to help understanding multiferroi'cs, mate- 
rials for spintronics applications, as well as temperature- 
dependent responses to magnetic fields, as needed for 
NMR experiments. 

Although we have focused on insulators at zero Kelvin, 
a generalization to metals at finite temperatures likely ex- 
ists, as in the case of the orbital magnetization (the n — 1 



term in our expansion), as outlined in Ref . [20I. Note, how- 
ever, that the density matrix idempotency relationship 
p = pp is not valid for such cases, so the density-matrix 
perturbation theory as developed in Ref. cannot be 
applied straightforwardly. 
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